Understanding Electron Transfer Reactions Using Constrained Density Functional Theory: Complications Due to Surface Interactions

The kinetic rates of electrochemical reactions depend on electrodes and molecules in question. In a flow battery, where the electrolyte molecules are charged and discharged on the electrodes, the efficiency of the electron transfer is of crucial importance for the performance of the device. The purpose of this work is to present a systematic atomic-level computational protocol for studying electron transfer between electrolyte and electrode. The computations are done by using constrained density functional theory (CDFT) to ensure that the electron is either on the electrode or in the electrolyte. The ab initio molecular dynamics (AIMD) is used to simulate the movement of the atoms. We use the Marcus theory to predict electron transfer rates and the combined CDFT-AIMD approach to compute the parameters for the Marcus theory where it is needed. We model the electrode with a single layer of graphene and methylviologen, 4,4′-dimethyldiquat, desalted basic red 5, 2-hydroxy-1,4-naphthaquinone, and 1,1-di(2-ethanol)-4,4-bipyridinium were selected for the electrolyte molecules. All of these molecules undergo consecutive electrochemical reactions with one electron being transferred at each stage. Because of significant electrode–molecule interactions, it is not possible to evaluate outer-sphere ET. This theoretical study contributes toward the development of a realistic-level prediction of electron transfer kinetics suitable for energy storage applications.

For a two-state system, the electron transfer rate constant (k ET ) is defined as where ν n and k n indicate nuclear vibration frequency and nuclear reorganization factor. In addition, k el is electronic transmission coefficient defines as It indicates the likelihood of an electron transfer taking place between a molecular acceptor and donor, or a transition (crossing) from the initial state I to the final state J. When the energy surfaces are parabolic and identical in shape, the Landau-Zener expression can be used to define P IJ . When the energy surfaces are parabolic and identical in shape can be defined by the Landau-Zener expression: 1,2 where ν el (S −1 ) is the electron hopping frequency in the intersection region at ta donor-acceptor separation r (Å), H IJ is the diabatic electronic coupling at separation distance r, λ is the reorganization energy, finally, h, k B , and T Planck's constant, Boltzmann constant, and temperature, respectively.
The electronic coupling element, H IJ , has an exponential distance dependence

S2
Taken together all the above equations, one can write where, ν 0 el , the value of ν el when r = r 0 , is defined as Combing Eqs. (S2) and (S6), one obtains . (S8) The analytical simplification of Eq. (S8) occurs based on two assumptions: This condition will be met for (i) outer-sphere ET (exp[−β(r − r 0 )] << 1) and (ii) nonadiabatic ET ( ν 0 el 2νn << 1). Using Taylor series expansion, we obtain Thus, Finally, one can write • On the other hand, if Similarly, the ET will be adiabatic or inner-sphere with: While for dBR5 and 2HNQ, graphene retains -2 electrons. . q X (e), q G (e), q total (e) denote total charges on the molecules, graphene, and whole system, respectively. In model (II), the missing charge is located on water molecules. The redox properties of the considered molecules can be described as a proton-coupled electron transfer (PCET) where the electron and proton transfer simultaneously or by a decoupled sequence of an electron transfer (ET; blue arrows in Figure S1) and a proton transfer (PT; green arrows in Figure S1) step. Readers interested in this topic are referred to the literature. 5,6 Combining these possible reaction paths results in a "scheme of squares" which is commonly used to evaluate and represent their redox properties ( Figure S1).
First, we aim to show that the protonation steps are unavoidable at a range of pH.
Therefore, we focus on the energetics of the decoupled electron and proton transfer reactions.
Both dBR5 and 2HNQ molecules are of the A type and thus, be fully reduced by taking two electrons and protons in the process: The energetics of the ET step, e.g. for A + e − ⇀ ↽ A − ., can be straightforwardly computed as:

S5
Where, G(A − ) and G(A) correspond to the Gibbs free energies of the reduced species and the oxidized states, respectively. G e − indicates the Gibbs free energy of an electron.
Similarly, the energetics of PT steps, e.g. for A + H + ⇀ ↽ AH + ., can be defined as Here, G(AH + ) corresponds to the Gibbs free energies of the protonated species, and G H + represents a proton's Gibbs free energy.
In order to compute the reaction of the PT and ET, we need to know the Gibbs free The pK a value is defined as: where R and T are the general gas constant and temperature, respectively. In addition, ∆G reaction indicates the difference between reactant and product Gibbs free energy. Using Eq. S18, the value of ∆G reaction and, consequently, G H + are determined (G H + = G reaction + G HCOOH − G HCOO − ). Therefore, the pK a value can be calculated for all the considered reactions.
We already know that cost of a proton used for redox reaction occurred in an implicit solvation model. Considering standard hydrogen electrode (SHE) as a reference electrode which is based on: the electron's Gibbs free energy is calculated as: Here, G H 2 (g) denotes the Gibbs free energy of H 2 in the gas phase. We calculate the Gibbs free energies of HCOOH and HCOO − compounds in the implicit water and H 2 in the gas phase. In order to determine G H + and G e − values, we use equations (S17), (S18), and (S19). Our calculation setup results in a value of −12.26 and −3.62 eV for G H + and G e − , respectively. To compute onset potential and pK a values, we first calculate the Gibbs free energies of reactant and product agents in the implicit water. Then, we compute reaction-free energy for ET and PT using equations (S15) and (S16), respectively. Finally, we define onset potential E 0 sol. , which corresponds to the required potential to involve only an electron in reaction, as: and pK a as Eq. S18. Additionally, the most stable tautomers were identified by computing the energies of the different structures where needed.

S7
faces Figure Figure S2 (a)). Contrarily, it occurs in the Normal region (see Figure S2 (b)). Single-minimum energy surface appears when the electronic S8 interaction is strong enough to be greater than the activation barrier-free energy (see Figure   S2 (c)). If not, a double-minimum energy profile is anticipated (see Figure S2 (d)).
5 Snapshot taken along the CDFT-AIMD Figure S3: Snapshot from CDFT-AIMD simulations when neutral OH-Vi molecule stands on top of neutral graphene.
It appears that water molecules are attempting to peal the OH-Vi out of the graphene after around 1.5×10 3 steps of CDFT-AIMD simulations for neutral graphene and the OH-Vi molecule. One snapshot is represented in Figure S3.
6 Finite-size effects and scheme of supercell extension A significant change in the Fermi-level energy in smaller supercells with respect to their charge state is shown in Figure S4 (a). As N , the number of carbon atoms in the supercell, increases, we observe that the variation decreases. Therefore, for the electrode to maintain a constant potential, the size must be infinite. However, we utilize a supercell of 15×9 unit cell containing 540 carbon atoms whose Fermi-level energy varies less than 130 meV between two successive charge states. In Figure S4 (b) and (c), the ∆A 0 IJ and λ values show systematic increases and decreases compared to the 7×4 graphene supercell, respectively. Those 15×9 graphene supercells are generated by embedding carbon atoms into the graphene from the 7×4 resulting CDFT-AIMD trajectories, as shown in Figure S5. While the atoms from step S9 Figure S4: (a) Fermi-level energy of graphene at different charges (q) versus the number of carbon atoms. Reference is made to neutral graphene. Parity plots of (b) reaction free energy (∆A 0 IJ ) and (c) reorganization free energy (λ) evaluated by 15x9 graphene (540 C atoms) and 7×4 graphene (112 C atoms). the latter calculations are frozen, the geometry of the added carbon atoms is optimized using DFT. Figure S5: Geometry of a 7×4 graphene (inner supercell) containing explicit water is increased to 15x9 by patching graphene unit cell.